home *** CD-ROM | disk | FTP | other *** search
/ Languguage OS 2 / Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO / gnu / glibc108.zip / glibc108 / sysdeps / ieee754 / ldexp.c < prev    next >
C/C++ Source or Header  |  1993-11-20  |  4KB  |  142 lines
  1. /* Copyright (C) 1992 Free Software Foundation, Inc.
  2. This file is part of the GNU C Library.
  3.  
  4. The GNU C Library is free software; you can redistribute it and/or
  5. modify it under the terms of the GNU Library General Public License as
  6. published by the Free Software Foundation; either version 2 of the
  7. License, or (at your option) any later version.
  8.  
  9. The GNU C Library is distributed in the hope that it will be useful,
  10. but WITHOUT ANY WARRANTY; without even the implied warranty of
  11. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12. Library General Public License for more details.
  13.  
  14. You should have received a copy of the GNU Library General Public
  15. License along with the GNU C Library; see the file COPYING.LIB.  If
  16. not, write to the Free Software Foundation, Inc., 675 Mass Ave,
  17. Cambridge, MA 02139, USA.  */
  18.  
  19. /*
  20.  * Copyright (c) 1985 Regents of the University of California.
  21.  * All rights reserved.
  22.  *
  23.  * Redistribution and use in source and binary forms are permitted provided
  24.  * that: (1) source distributions retain this entire copyright notice and
  25.  * comment, and (2) distributions including binaries display the following
  26.  * acknowledgement:  ``This product includes software developed by the
  27.  * University of California, Berkeley and its contributors'' in the
  28.  * documentation or other materials provided with the distribution and in
  29.  * all advertising materials mentioning features or use of this software.
  30.  * Neither the name of the University nor the names of its contributors may
  31.  * be used to endorse or promote products derived from this software without
  32.  * specific prior written permission.
  33.  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
  34.  * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
  35.  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
  36.  */
  37.  
  38. #include <ansidecl.h>
  39. #include <math.h>
  40. #include <float.h>
  41. #include "ieee754.h"
  42.  
  43. __CONSTVALUE double
  44. DEFUN(ldexp, (x, exp),
  45.       double x AND int exp)
  46. {
  47.   union ieee754_double u;
  48.   unsigned int exponent;
  49.  
  50.   u.d = x;
  51. #define    x u.d
  52.  
  53.   exponent = u.ieee.exponent;
  54.  
  55.   /* The order of the tests is carefully chosen to handle
  56.      the usual case first, with no branches taken.  */
  57.  
  58.   if (exponent != 0)
  59.     {
  60.       /* X is nonzero and not denormalized.  */
  61.  
  62.       if (exponent <= DBL_MAX_EXP - DBL_MIN_EXP + 1)
  63.       {
  64.       /* X is finite.  When EXP < 0, overflow is actually underflow.  */
  65.  
  66.       exponent += exp;
  67.  
  68.       if (exponent != 0)
  69.         {
  70.           if (exponent <= DBL_MAX_EXP - DBL_MIN_EXP + 1)
  71.         {
  72.           /* In range.  */
  73.           u.ieee.exponent = exponent;
  74.           return x;
  75.         }
  76.  
  77.           if (exp >= 0)
  78.           overflow:
  79.         {
  80.           CONST int negative = u.ieee.negative;
  81.           u.d = HUGE_VAL;
  82.           u.ieee.negative = negative;
  83.           errno = ERANGE;
  84.           return u.d;
  85.         }
  86.  
  87.           if (exponent <= - (unsigned int) (DBL_MANT_DIG + 1))
  88.         {
  89.           /* Underflow.  */
  90.           CONST int negative = u.ieee.negative;
  91.           u.d = 0.0;
  92.           u.ieee.negative = negative;
  93.           errno = ERANGE;
  94.           return u.d;
  95.         }
  96.         }
  97.  
  98.       /* Gradual underflow.  */
  99.       u.ieee.exponent = 1;
  100.       u.d *= ldexp (1.0, (int) exponent - 1);
  101.       if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
  102.         /* Underflow.  */
  103.         errno = ERANGE;
  104.       return u.d;
  105.       }
  106.  
  107.       /* X is +-infinity or NaN.  */
  108.       if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
  109.       {
  110.       /* X is +-infinity.  */
  111.       if (exp >= 0)
  112.         goto overflow;
  113.       else
  114.         {
  115.           /* (infinity * number < 1).  With infinite precision,
  116.          (infinity / finite) would be infinity, but otherwise it's
  117.          safest to regard (infinity / 2) as indeterminate.  The
  118.          infinity might be (2 * finite).  */
  119.           CONST int negative = u.ieee.negative;
  120.           u.d = NAN;
  121.           u.ieee.negative = negative;
  122.           errno = EDOM;
  123.           return u.d;
  124.         }
  125.     }
  126.  
  127.       /* X is NaN.  */
  128.       errno = EDOM;
  129.       return u.d;
  130.     }
  131.  
  132.   /* X is zero or denormalized.  */
  133.   if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
  134.     /* X is +-0.0. */
  135.     return x;
  136.  
  137.   /* X is denormalized.
  138.      Multiplying by 2 ** DBL_MANT_DIG normalizes it;
  139.      we then subtract the DBL_MANT_DIG we added to the exponent.  */
  140.   return ldexp (x * ldexp (1.0, DBL_MANT_DIG), exp - DBL_MANT_DIG);
  141. }
  142.